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Abstract 

We review some recent developments of Grad's approach to solving the Boltzmann 
equation and creating reduced description. The method of invariant manifold is 
put forward as a unified principle to establish corrections to Grad's equations. A 
consistent derivation of regularized Grad's equations in the framework the method 
of invariant manifold is given. A new class of kinetic models to lift the finite-moment 
description to a kinetic theory in the whole space is established. Relations of Grad's 
approach to modern mesoscopic integrators such as the entropic lattice Boltzmann 
method are also discussed. 

1 Introduction 

There has been a long-standing quest for improving on the Grad 13-moment approxima- 
tion [34J. In particular, such an improvement is needed to study the interplay between 
hydrodynamics and kinetics in the domain of moderate Knudsen numbers, in particular, 
simulations of flows at a micrometer scale in so-called micro-electro-mechanical systems 
(MEMS) ^T]. The recent renewed interest to this topic is consistent with the current 
trend in computational fluid mechanics to use minimal kinetic models instead of more 
traditional numerical schemes for hydrodynamic equations. 

Let us recall the famous Grad's 10-moment and 13-moment approximations for the 
distribution function: 
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Here, as usual, f^^^ is the local Maxwellian, c = v^^{v — u) is the "peculiar" velocity, u 
is the local flow velocity, vt = \Jm/2k-QT is thermal velocity, p is scalar pressure, cr is the 
nonequilibrium stress tensor, and q is the heat flux. 
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Technically, in Grad's original approach, parametric families (P) and Q were intro- 
duced as truncated Hermite polynomial expansions of the distribution function around 
local Maxwellians. However, it is much more attractive to view Grad's distributions as 
parametrically specified sub-manifolds ( "surfaces" ) in the larger space of distribution func- 
tions. Grad's method has given start to a host of new methods focused around the hard 
question of nonequilibrium statistical mechanics: How to effectively reduce the microscopic 
to a macroscopic description? This review is devoted to some selected instances of this 
question. 



2 Grad's method and beyond 

"Grad's legacy" (where and how to go beyond the 13-moment approximation?) was inter- 
preted and extended in different ways by many authors. Let us mention those which are 
most relevant to the present discussion. 



Quasi-equilibrium approximation 

Quasi- equilibrium approximation (or maximum entropy approximation) in the application 
to the Boltzmann equation was established in the sixties by several authors, in partic- 
ular, by Kogan jlHl and Lewis |1HI, though we note that it was already mentioned by 
Grad himself, and also by Koga (cf. |1H])- A detailed discussion geometrical aspects of 
quasi-equilibria was given in jHS]. The construction is based on solving the conditional 
maximization problem: For the concave functional S = —k^ J fin fdv (local entropy den- 
sity) and for given distinguished linear functionals M{f), find 

S max, M{f) = M. (3) 
The solution in terms of Lagrange multipliers (dual variables) A is written as 

/ = exp|^AfcZ}yM,|. (4) 
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If Mk{f) = J mk{v)fdv, then we have 



/ = exp <^ ^ AkTUk 



(5) 
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If M = Mq = J{l,v,v'^}fdv, the parametric set (jlj) coincides with the set of local 
Maxwellians. If M = /{I, V, vv}fdv this is the 10-moment quasi-equilibrium approxima- 
tion, whose expansion to linear order in cr/p coincides with Grad's 10-moment approxi- 
mation Though almost pedantic, some attention is required when proceeding to the 
13-moment case: Functions Q are defined in terms of dual variables. What is not always 
well defined is the moment chart (or moment parameterization) of these sets, or, A(M), and 
a regularization of divergent integrals is required. This can be done either by restricting 
the velocity integration domain to a large ball v < y/E, where E is the total kinetic energy 
of the gas in a container ^E], or by introducing a higher-order even velocity polynomial 
(at a price of an extra variable). In the first case of regularization it is possible to use the 
smallness of q/pvx to expand the regularized distribution and send the radius of the ball 
to infinity to end up with the 13-moment Grad's approximation (j21). 

A particularly useful version of entropic methods was introduced in [39|| for chemically 
reacting gas mixtures, and discussed in some detail for a single-component gas in [THl |21] 
[triangle entropy method). It can be viewed either as a stepwise realization of the basic 
maximization problem, or (better) as a self-consistent recipe. Let us split the totality of 
distinguished macroscopic variables M into M' and M" , where M' are linear functionals 
for which we can solve explicitly the problem Q ("easy variables"), and where M" are 
"difficult variables". The difficult variables may be even nonlinear in the distribution 
function, for example scattering rates (see below), so that even the statement of the problem 
of the entropy maximization can cause difficulties. The triangle entropy method allows to 
construct quasi-equilibria even in these cases. Let us denote f{M') the quasi-equilibrium 
corresponding to the easy variables. Then the triangle quasi-equilibrium for the difficult 
variables is found as follows: expand the entropy functional up to quadratic order around 
f{M'), and find maximum under the conditions that (i) Easy variables M' are not changed, 
and (ii) Difficult variables are fixed to linear order. That is. 



Maximization here is with respect to 6f, nonlinear parametric dependence on M' is not 
varied. The nice property is that © is always solvable in closed form, the resulting triangle 
quasi-equilibrium. 




(6) 



/(M', 5M") = f{M') + 6f{M', 6M"), 



(7) 
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depends linearly on 5M" and nonlinearly on M' (so the overall dependence is quasi-linear) . 
If M' are the five hydrodynamic fields, and if M" are J vvfdv or J {vv, vv'^}fdv (both easy 
and difficult variables are hnear in this example), then Q are Grad's 10- and 13- moment 
distributions, respectively |18j . 

The advantage of the quasi-equilibrium approximations is that they are equipped natu- 
rally with the thermodynamic parameterization [201 ■ The structure of the thermodynamic 
parameterization assumes specification of the projector P onto the tangent bundle of the 
quasi-equihbrium manifold. For quasi-equihbria, this is (we stick to the case (jSJ for sim- 
plicity): 

k 

The purpose of projector P is to define dynamics along the manifold. Namely, if we 
write the Boltzmann equation, 

Dtf = J{f) = -{v-u)-Vf + Qif, /), (9) 

where Q{f, /) is the Boltzmann collision integral, and Dt = dt + u-'V is the time derivative 
in the co-moving reference system, then the vector field attached to each state on the quasi- 
equilibrium manifold (or the microscopic time derivative) is: 

Dr-/(M) = J(/(M)). (10) 

Here, we simply evaluate the action of the operator in the right hand side of the Boltzmann 
equation (jH)) on the quasi-equilibrium distributions. On the other hand, under the action 
of the projector P (jSI), vectors J{f{M)) yield the vector field on the tangent bundle of the 
quasi-equilibrium manifold, or the macroscopic time derivative: 

A™/(M) = PJ{f{M)). (11) 

The latter can be viewed as a short-hand writing of Grad's equations, which follow from 
flll|) upon multiplication with m^, and integration: 

dtMk + u ■ VM, = j mkPJ{f{M))dv. (12) 

One can ask, what is the use of the microscopic time derivative fjlOp when only its 
projected piece, PJ(/(M)), contributes finally to Grad's moment equations (fT^? The 
answer is that the comparison of the vectors J(/(M)) and PJ(/(M)) measures how good 
the closure (fT^ really is. The difference between J(/(M)) and PJ(/(M)) is of such a 
great importance that it deserves a specific name. The defect of invariance (of the quasi- 
equilibrium approximation) is, 

A{M) = J{f{M))-PJ{f{M)). (13) 
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A moment representation of the defect is also useful: If mi, ... , m„ are the distinguished 
moments, and . . . are the higher-order moments, then the velocity-dependent func- 

tion is equivalent to the infinite sequence Aj(M) by taking moments of ()13|l : 



Levermore proved hyperbolicity of maximum entropy approximations. Dual parame- 
terization of quasi-equilibrium manifolds prove to be an advantage in numerical realizations 
(so-called Legendre integrators, see e. g. [^EnilHl in the context of polymer dynamics). 

The quasi-equilibrium approximations reveal most clearly the time hierarchy assump- 
tion behind the Grad's approach jlHI : In the fast relaxation, the entropy grows according to 
Boltzmann equation until maximum of entropy is reached on the quasi-equilibrium states. 
After that the slow evolution takes place along the manifold of the quasi-equilibrium states. 
Putting this assumption on trial and, if needed, improving on it is the key in seeking cor- 
rections to Grad's approximation. The trial is the deviations away from zero of the defect 
of invariance (fT^ . 

In this section, we reviewed the basic structure of Grad's theory, and indicated that 
a way beyond a given moment approximation should take into account the defect of its 
invariance. Before proceeding along this line, let us discuss two other routes, which can be 
indicated as "increase the number of variables" and "take other variables" . 

Many moments approximations 

With a given moment approximation at hand, and without asking the question, "How good 
is this approximation?" there is only one option to try to improve on it - to extend the 
list of distinguished variables, and to construct another approximation. This viewpoint 
dominated earlier studies on moment approximations, and was followed by many authors, 
in particular, by Miiller and Ruggeri ^] and their collaborators, mostly in the quasi- 
linear form, using orthogonal functions developments, and for the Bhatnagar-Gross-Krook 
model collision integral with phenomenological dependence of the relaxation parameter. 
Convergence to Boltzmann equation is a difficult question; in fact, it is not expected that 
Grad's distributions converge to solutions of the Boltzmann equation, even pointwise J2] • 
On the other hand, as numerical results show [HI], the weak convergence (convergence of 
the moments) can be expected in the linear case. However, without at least evaluating 
the defect, Eqs. fll3|) or ()14|) . the uncontroUability of approximations with any number of 
moments remains. 

Scattering rates as independent variables 

Remarkably, Grad's approximations with moments as slow variables (0) and (j21) (or any 
other with more moments) do not contain any molecular information. This shortcoming is 
inherited from the use of simple sets of orthogonal functions in the original Grad method. 




0, i = 1, . . . ,n 



J m,{J{f{M)) - PJ{f{M)))dv, z = n + 1, . . . 



(14) 
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However, if one thinks of using less variables to capture more physics, then other (non- 
moment) quantities can be tried. In particular, for the pure gas, interesting variables are 
scattering rates of moments, for example. 



o-^°ii = j mvvQifJ)dv, (15) 

is the scattering rate of the stress tensor, Q is the Boltzmann collision integral. The 
variable cr™^\ unlike cr, contains information about the molecular interaction. Using the 
triangle entropy method with the five hydrodynamic variables as easy, and cr^°'' as difficult 
(this is indeed the case because cr'^°^^ is nonlinear in /), we can construct a "scattering" 
counterpart of the 10- moment approximation (^: 

where ^o{T) is the first Sonine polynomial approximation to the viscosity coefficient, and 
the difference from Grad's 10-moment approximation is in the dimensionless function R. 
The function R depends on the particle's interaction, and R is constant only for Maxwell 
molecules in which case the present approximation is equivalent to Grad's approximation 
(up to renaming the variables). For hard spheres j.24j . 



R 



j\-'''''{i-f)(^(i-e) + 2)dt. (17) 

The case when the scattering rate of the heat flux is included (the counterpart of the 
13-moment approximation), and also a mixed version (moments and scattering rates both 
as distinguished variables), were studied in Ref. [TH] . 

Eventually, the triangle entropy method makes it possible to handle "smart" variables 
which may be more appropriate to the physics of the problem at hand rather than plain 
moments. The latter assumes a certain degree of a physical intuition; furthermore, the 
uncontrollability of the resulting quasi-linear quasi-equilibria remains an issue. 



Method of invariant manifold 

The general method to derive dynamic corrections on top of successful initial approxima- 
tions like Grad's was developed by the authors [201 UHl 1211 1101 12Z|- The essence of 
the method of invariant manifold (MIM) is (i) to write the invariance condition in the 
differential form (the microscopic time derivative on the manifold equals the macroscopic 
time derivative), and (ii) to solve this equation by iterations. The choice of the initial 
approximation is an important problem. Often, it is convenient to start from the quasi- 
equilibrium manifold (this will be our choice below in this paper). However, the choice of 
the initial manifold in MIM is not restricted to quasi-equilibria. The typical example gives 
us the famous Tamm-Mott-Smith approximation for the strong shock wave (see discussion 
in [201123123 )• Strictly speaking, the method of invariant manifold can be applied in order 
to reflne any initial approximation compliant with some transversality conditions. 
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Correction to local Maxwell manifold 



In order to explain the two steps of the MIM, we shall consider first, for the purpose of 
illustration, the invariance correction to the local Maxwell approximation to all orders in 
Knudsen number. 

The main idea is to pose the problem of a finding a correction to Euler's hydrodynamics 
in such a way that Knudsen number expansions do not appear as the necessary element of 
the analysis. This will be possible by using Newton method instead of Taylor expansions 
to get such correction. 

The starting point is the manifold of local Maxwell distribution functions (LM) /o(n, u, T; v) 
where v is particle's velocity, and n, u, and T are local number density, average velocity, 
and temperature. As appropriate to our approach, we first check the invariance of the LM 
manifold. Projector (jSJ on the LM manifold is: 



fo 

n 



Jdc + 2c • y cJdc + 3 (^c^ - 2 J j [f ~ 2) "^^^ ' ^'^^^ 

Computing the microscopic time derivative ()10|) on the LM states, projecting it with Pq 
(fTS|) to get the macroscopic time derivative (the time derivative due to Euler's equations) 
(fTT| . and subtracting the second out of the first, we evaluate the invariance defect of the 
LM manifold 



A(/o) = Afo) - PoJ{h) = -fo 



2Vu : I cc- -Ic" 



VT 
+ vt— ■ c 



(19) 



The defect is not equal to zero as long as there average velocity and the temperature vary 
in space, as expected. Note that the defect is neither small or large by itself. 
To find the correction to the LM manifold, we write the invariance condition, 

A(/) = J{f) - PJif) = 0, (20) 
and consider it as an equation to be solved with the initial approximation /q for the 
manifold and Pq for the projector since both are unknown in the (j20|) . This might seem too 
much to require, however, the well-posedness of the problem is restored once the additional 
requirement that the manifold we are looking for should be the manifold of slow motions 
is invoked (see Ref. [23J for details). Here we will consider the first iteration. 

Upon substitution of Pq in place of P, and of /i = /o + Sf in place of / in equation 
(1201), after the linearization in 6f, we get 



Lf.Sf + (Po -l){v-u)- V6f + A(/o 



0, (21) 

where L/^ is the linearized collision integral (linearization in the local Maxwell state, and 
we keep indicating the linearization point for reasons to be seen later). Equation (PT|) has 
to be solved subject to the condition. 



PoSf = 0. 



f22) 
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The linear equation of the first iteration ()2H) is the most important object in our theory. 
Indeed, it does not contain at all the smallness parameter, and, in fact, was obtained 
without assumption of a small Knudsen number. If, however, the Knudsen number is 
introduced into equation by the usual rescaling, 

Lfo -Lfo, (23) 
then the first-order in e solution, 6fi ~ e6f[^^ is found from the integral equation, 

Lf,6fi'^ = -A(/o), (24) 

which has to be solved subject to the condition Pq^/i^''- K is readily checked that equation 
()24|) is just the equation of the first approximation of the Chapman- Enskog method, which 
is thus recovered as the special case of MIM in the collision-dominated limit. 

The invariance equation for the first correction (|2ip is more complicated than the 
Chapman-Enskog equation because it also contains the spatial derivatives (though 
it is much less complicated than the linearized Boltzmann equation because the there is 
no time dependence in the equation (jUJ). Methods to treat equations of the type 
has been also developed in [101 EHl I2Z| and worked out for many kinetic systems (not 
only for the Boltzmann equation). Here we consider, for the purpose of illustration, the 
small deviations around the global equilibrium in 3D. We shall treat equation 1)211) in 
such a way that the Knudsen number will appear explicitly only at the latest stages of the 
computations. 

We denote as Fq the global equilibrium with the equilibrium values of the hydrodynamic 
quantities, uq, Uq = 0, and Tq. Deviations are 6n, 6u, and 6T. We also introduce reduced 
deviations. 

An = Sn/no, An = 5u/v^, AT = ST/Tq, 

where is thermal velocity in the equilibrium. 
We seek the invariance correction, 

h=Fo{l + ipo + ipi), (25) 

where 

^0 = An + 2Au ■c + At(^c^-^. (26) 

comes from the linearization of the local Maxwellian around Fq, and where (fi is the 
unknown function to be found from equation invariance equation ()2H). In order to find (pi, 
we apply a Galerkin approximation in order to achieve a finite-dimensional approximation 
of the linear collision operator j^, which amounts to setting 



^^ = A{x) .c[c'-^^+ B{x) : (cc - he' ) . (27) 
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Our goal is to derive functions A and B from a linearized version of equation ()2H). 
Knowing A and B, we get the following expressions for shear stress tensor cr and heat flux 
vector q: 

a = poB, q = ^Pov'tA, (28) 

where po is the equilibrium pressure. 

Linearizing equation (pT|l in Fq, substituting ipi ()27|1 . and switching to the Fourier- 
transformed in space variables, we derive the set of linear algebraic equations for the 
Fourier image of the functions A and B (which we denote as a^, and bk, respectively): 



-iv^Tkn; (29) 

where i = y/—T, fc is a wave vector, /iq is the first Sonine polynomial approximation of the 
shear viscosity coefficient, Tk and 7^ are Fourier images AT, and Au, respectively, and the 
over-bar denotes a symmetric traceless dyad. 

Introducing dimensionless the reduced wave vector, 

Po 

solution to equation may be written: 

3[(5/3) + {1/2)k^] + 3[(5/3) + (l/2)/2][5 + 2/€2] ~ 2[5 + 2/€2]' 

_ 5[/^(7fc-/^)+7fc/'(5 + 2/t^)] 
2[5 + 2/t2] 3[5 + 2/€2][(5/3) + (l/2)fi:2] ' ^ > 

With the Fourier- image of the fluxes ()28p. 

5 

o-fc =Pobk, q = -poVTdk, 

which have to be used to close the Fourier-transformed linear balance equations, func- 
tions ()30|) concludes our computation of the dynamic correction to the linearized local 
Maxwellian. Note that due to the non-polynomial in k contributions, the resulting lin- 
ear hydrodynamics is highly nonlocal. This is, of course, not surprising because no small 
Knudsen number expansions truncated to some order ever appeared. 

Let us briefly consider the new hydrodynamic equations specializing to the one-dimensional 
case. Taking the z-axis for the propagation direction, and denoting kz as k, 7 as 7^, we 
obtain in ()30|) the Fourier images of a = a^; and b = bzz (full notation are restored here): 



5po 
3/io 



ttfc + iv^bk ■ k 



— 6fc + iv^ka^ 
Mo 



bk = 
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1 + VgK)^^^ 

1 + |pn • 



5P0 ^ol^T 

These expressions close the hnearized balance equations 



■\dtPk + ik'jk = 0, (32) 



2 

3 5 



In order to restore the Knudsen number in (jSU, we introduce ^m.f.p. = v^fJ'o/Po (the 
quantity /m.f.p. is of the order of the mean free path), and we also introduce a hydrodynamic 
scale Ih, so that k = n/lh, where k is dimensionless. With this, we obtain in the equation 



b 



l + le^K^ ' ^^^^ 
4 • I 2 2 



K 

I -I- ^ 



1 + h^K^ 



where e = /c/^/i is the Knudsen number. In the limit e — > 0, equation reduces to the 
familiar Navier-Stokes-Fourier expressions: 

4 

(^zz = --fJ'odzSuz, qz = -XodJT 

where Aq = 15kBfio/4:m is the first Sonine polynomial approximation of the thermal con- 
ductivity. 

In order to find out a result of non-polynomial behavior ()33|1 . it is most informative to 
calculate a dispersion relation for planar waves. It is worthwhile introducing dimensionless 
frequency A = ulh/v^, where u; is a complex variable of a wave ~ exp(ci;t + ikz) (Reu is 
the damping rate, and Imu; is the circular frequency). Making use of ()32|) and (jH^ . writing 
e = 1, we obtain the following dispersion relation \{k,): 

12(1 + "^n^X^ + 23k\1 + -K^)X'' + 2k\5 + + + + -k^) = 0. (34) 

5 5 5 2 5 
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Figure 1: Attenuation rate of sound waves. Dotts: Burnett approximation. Bobylev's 
instability occurs when the curve intersects the horizontal axis. Solid: First iteration of 
the Newton method on the invariance equation. 

Fig. ^ presents a dependence ReA(K^) for acoustic waves obtained from (jMj) and for 
the Burnett approximation ^|. The violation in the latter occurs when the curve crosses 
the horizontal axis. In contrast to the Burnett approximation ^2], the acoustic spectrum 
(IH4|1 is stable for all k. Moreover, ReA demonstrates a finite limit, as k ^ oo (so-called 
"Rosenau saturation" [ST]). 

The example considered demonstrates how to apply the method of invariant manifold 
in the simplest case of the initial manifold. Let us now switch back to another case of the 
initial manifold, the Grad's 13-moment approximation. 

Invariance correction to 13-moment manifold 

As we've said before, MIM is able to address the invariance correction, in principle, to 
any interesting initial approximation, so it may be not surprising that the next candidate 
after the local Maxwell manifold are manifolds of Grad's distributions. The problem of 
finding the invariance correction to the moment approximations was first addressed in Ref. 
IHj. Without repeating computations of the Ref. [Sj, our objective here is to explain 
why the correction to the Grad's manifolds is a distinguished case. Let us start with 
the quasi-equilibrium of a generic form, and compute the first iteration of the invariance 
equation, 

(1 - P)Lf(M)6f + {P-l){v-u)- V6f + A(M) = 0, (35) 

where P is the quasi-equilibrium projector (jSJ, A(M) is the invariance defect of the 
quasi-equilibrium approximation, and Lf(^M) is the collision integral, linearized in the quasi- 
equilibrium. Notice that the latter object is not well studied in the classical theory of the 
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Boltzmann equation, where most of the reduction problems are using the colhsion integral, 
linearized in the local equilibrium. A simplification of equation ()35p happens when we look 
for Grad's quasi-linear approximations, then, to the linear order accuracy in the higher- 
order moments, we can replace, 

Lf{M) Lfg, (36) 

in the equation ()35|). where Lj^ is the usual linearized collision integral in the local 
Maxwell state. 

Let us proceed further with evaluation of the other terms in the equation (j^ . The 
projector, corresponding to the 13-moment Grad's approximation, reads, 

Pi3 = Po + ni3, (37) 
where Pq is the local-equilibrium projector (fTH|l . and where Ilia acts as follows: 

Ui3J=^!^Y : j YJdv + Z- J ZJdv"^ , (38) 

Here Y = -\/2cc, and Z = -^c (c^ — |) , are peculiar velocity polynomials forming the 
13-moment approximation. 

Computing the defect of invariance of the 13 moment approximation to the linear order 
in the if 13, we see, that there are two contribution, local (containing the linearized collision 
integral), and the nonlocal (containing the free flight operator). 



+ (39) 
{l-Urs)[LfJoipn], 
(l-P,,)[-{v-u)-\/foil + ipi3)]. 

Before proceeding any further, we shall discuss the physical significance of the defect 
(jnni) because it is the first instance where classical methods, like the Chapman-Enskog 
method, become inapplicable. 

The nonzero defect of invariance of any manifold reveals the following: The solution to 
the Boltzmann equation with the initial condition on the manifold leaves this manifold at 
t > 0. The two parts of the defect correspond to two different mechanisms responsible for 
this to happen. The local defect is not equal to zero whenever the polynomials Y and Z, 
forming the Grad manifold, are not eigenvectors of the linearized collision integral. This 
is distinct from the dynamic noninvariance of local Maxwellians, where in the latter case 
the local defect is equal to zero whatever the collision model is chosen. For the Grad 
approximation, A'^°3^ = in only two (commonly known) cases, i. e. for Maxwell molecules 
and for the Bhatnagar-Gross-Krook model (BGK). 

It is important to recognize that, whenever the local defect is not equal to zero, the 
initial manifold has to be first corrected locally, in order to bring it closer to the slowest 
eigenspace of the collision operator, and before any nonlocal corrections due to A°3°'^ are 



^13 
A loc 
^13 
A nloc 
^13 
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addressed. It has been demonstrated in 03] that the first local correction to the 13- 
moment approximation results in Grad's equations corrected in the sense that the transport 
coefficients become the exact Chapman-Enskog coefficients (not the ffist Sonine polynomial 
approximation, as in the original Grad equations). Whether or not the local correction 
is spectacular in the context of the single-component gas with traditional collision models 
like hard spheres (the ffist Sonine polynomial approximation is "good" already...), the 
clear distinction between local and nonlocal corrections is crucial. For example, with this 
distinction it is possible to extend the method of invariant manifolds to driven systems (see 
the derivation of the Oldroyd constitutive equations from polymer kinetic theory 62]). 

Now, let us turn our attention to the nonlocal piece of the defect of invariance. It 
can be demonstrated that A'f^'^ contains no terms with gradients of neither of the five 
hydrodynamic fields, the only gradient contributing to A^^'^ are of the stress tensor and of 
the heat flux. In the linear approximation near the global equilibrium Fq |44j . 

^n" = -v^rFoi^ilkrsdkCTrs + Tl2\iAql + Ilsdkqk) (40) 
where di = d/dxi, and 11 are velocity polynomials: 

^l\krs = Ck [CrCs - (1/3)(5^^C^] - (2/5)4^0^0^, 

n2Kfc = (4/5) [c^ - (7/2)] {ack - (l/3)5ifcc2] , 
ns = (4/5) [c^ - (5/2)] [c^ - (3/2)] - c^. 

The absence of the gradients of the hydrodynamic fields in the nonlocal defect reveals 
some important information: The invariance correction to the 13-moment approximation 
differs from the higher-order corrections to the hydrodynamic equations. For example, since 
the linearized hydrodynamic equations following from the 13-moment Grad (uncorrected) 
equations to second order in Knudsen number are the Burnett equations for Maxwell 
molecules (see, e. g., comparison of corresponding dispersion relations in |1S]), the same is 
true also for the corrected Grad equations, as explicitly verified by Struchtrup and Torillhon 

After this qualitative analysis of the defect of the invariance of Grad's approximation, 
let us finish setting up the invariance equation of the ffist iteration formally. With the 
replacement (jH^ in the equation and using Pq-^/o = 0; '^^ have, 

(1 - ni3)L^o5/ + (Pi3 -\){v-u)- V6f + Al- + A?^- = 0. (41) 

In principle, this equation can be studied in the same spirit as the equation of the ffist 
iteration to local Maxwellians, that is, without introducing small parameters. However, it 
is much more instructive to consider the collision-dominant case, introducing the scaling 

(ESI), 

(1 - U,s)hj,6f + (Pi3 -l){v-u)- V6f + ^Al°3^ + Afr = 0. (42) 
The correction 5f to ffist two orders, 5f ~ 6f^^^ + e6f^^^ is found from equations: 
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(1 - UnWf^'^ = (43) 
(1 - ni3)%//(^) = (44) 

These equations have to be solved subject to the additional conditions, Pi^Sf^^^ = 0, and 
Pis^f^^^ = 0, respectively. 

The first equation (jlH]) is responsible for the local correction, as expected. It's signif- 
icance was discussed above. For the following, we assume no local correction is needed, 
that is, either we assume BGK or Maxwell molecules (then 5f^^^ = rigorously), or that 
the Grad's approximation is a reasonably good approximation for the eigenvectors of Lf^. 

5/(1) = ^A5^^°^ (45) 

As it was emphasized in Ref. |13], the invariance correction to the 13-moment Grad's 
approximation is related to the 13-moment Grad equations entirely in the same way as the 
Navier-Stokes equations are related to the Euler equations. Roughly speaking, it uses the 
same amount of the Boltzmann collision integral as the classical first-order equation of 
the Chapman- Enskog method. For that reason, this is the distinguished case among other 
possible applications of MIM to improve on moment approximations. 



Strongly nonlinear invariance corrections 

As we have seen it in the previous section, the invariance correction to the quasi-linear 
quasi-equilibria (Grad's moment approximations) is distinguished by the fact that we can 
compute it with the usual linearized collision integral in the local equilibrium. Then the 
resulting linear integral equations have the same structure as in the classical Chapman- 
Enskog method, Lf^^p = A (albeit with a different right hand side A). The operator Lf^ is 
self-adjoint in the scalar product generated by the second differential of the entropy in the 
local equilibrium, and thus it is simple to solve. This is not the case when the manifolds we 
want to correct contain pieces well beyond the vicinity of the local equilibrium, for example, 
general quasi-equilibria. In these cases, the linearized collision operator Lj^m) is not self- 
adjoint anymore. In such cases, it was suggested to use the symmetric linearization in 
order to establish dynamic corrections in highly nonequilibrium situations. The symmetric 
linearization of the Boltzmann collision integral in the state / has the form. 



I f'fi + ffl 



/' f[ fl f 



dv[dv'dvi. (46) 



Note that Ly'™ ~^ -^/o if the state / tends to the local Maxwellian /o (the consequence 
of the detail balance, f f[ = f fi in the local equilibrium). Operator Lj^ enjoys the 
familiar properties of the usual linearized collision integral. Let us introduce notation for 
the entropic scalar product in the state /: For two distribution functions gi and g2, 

{gi\g2)f = J -y^^ — (47) 
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The following three properties of the operator L^J™ are immediate consequence of the 
definitions (gSl) and (gZj): 

(i) (fi'i|^7"'lfi'2)/ = (5-2 1 1 5-1 )/ (symmetry); 

(ii) {g\VJ"^\g) f < (local entropy production inequality); 

(iii) LY'^g = if g/f G Lin{l, v,f^} (conservation laws). 

Using the symmetric linearization, we see the equation for the invariance correction for 
a general quasi-equilibrium f{M) becomes 

(1 - P)Lj^l)5f +{P-l)(v-u). V5f + A(M) = 0. (48) 

The only difference with the equation ()35|1 is in the replacement of the linearized collision 
integral -Z^/(m) with the symmetric linearized collision integral L^^^y This difference is 
crucial though: When using symmetric operator we get back all the familiar tools for 
solving integral equations (Fredholm alternative in the collision-dominated case ^Hj, the 
parametrix expansion without such domination [SHI, and like). All this is impossible with 
the plain linearized operator I>/(m), for example, even the null-space of -Z^/(Af) is not known 
in general. 

The symmetric iteration was tested in the case of finite-dimensional kinetic systems of 
chemical kinetics |32[ |2Hj (see, in particular, its recent application to construction of grid 
representations of invariant manifolds jSI])- Results of convergence of symmetric iterations 
to slow nonlinear manifolds are quite encouraging. At the time of this writing, symmetric 
iteration remains almost entirely unexplored for the Boltzmann equation. 

Invariance principle in the moment representation 

The aforementioned computations of various quasi-equilibria, moment approximations and 
their invariance corrections were all done in the setting of the kinetic Boltzmann equation 
and distribution functions. Invariance corrections can also be studied in a simplified setting: 
Consider a closed system for n = k + m moments and reduce it to a closed system for k of 
them. Such a simplification (with respect to the full kinetic theory) makes sense especially 
if one wants to get a basic qualitative understanding about the form of reduced description 
in terms of k moments. 

This problem was studied to some very detailed extend, and well beyond the usual first- 
order Knudsen number corrections, for the case of hydrodynamics from 10- and 13-moment 
Grad equations beginning with the paper ^H] on a partial summation of the Chapman- 
Enskog expansion to all the orders in Knudsen number, and on the exact summation of the 
expansion [22] • Some of these studies were recently summarized in jlH! with the emphasis 
of the iteration method for solving the invariance equation, and the interested reader is 
directed to that paper. In spite of a seemingly drastic simplification with respect to the 
"true" kinetic theory, results are sometimes surprisingly robust. For example, the leading 
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invariance correction to the nonlinear longitudinal viscosity in so-called homoenergetic flow 
found from the 10-moment equations |H] is exactly the one found independently from the 
exact solution to the model Boltzmann equation jHSI- New interesting results on summation 
of the Chapman-Enskog expansion for the Boltzmann equation were obtained recently by 
Slemrod 

The framework of a larger moment system used to obtain the invariance correction to a 
smaller moment system appears in the recent work of Struchtrup and Torrilhon [32] ■ The 
difference with the derivation from the Boltzmann equation is basically the absence of the 
local correction. The study [59j demonstrated a set of advantages of these equations above 
the Grad's system, most importantly, the improved shock wave structure. 



3 Quasi-equilibrium kinetic models 

The invariance corrections explore more of the phase space than initially assumed by 
making a Grad approximation with a given number of variables. By measuring the defect of 
invariance, we realize in which direction the quasi-equilibrium manifold should be improved 
in order to take into account fast motion towards it. There is another useful way to explore 
fast motions: To lift the dynamics from the manifold to a dynamics in the full space by 
means of a kinetic model. 

We recall that lifting the Euler dynamics which takes place on the local Maxwell mani- 
fold to a kinetics in the whole phase space is done by the very useful Bhatnagar-Gross-Krook 
model (BGK), 

dJ + vVf = --{f-fo{f)), (49) 
r 

where r > is the relaxation time, and /o(/) is a map / — > /o established by local 
conservation laws: 

l{l,v,v'}if-fo{mv = 0. (50) 
The right hand side of Eq. , 

Qbgk = --(/-/o(/)), (51) 
r 

is called the BGK collision integral. Proof of the iZ-theorem for the BGK kinetic equation 
does not rely anymore on the microscopic reversibility (as in the Boltzmann case), instead, 
it is a simple consequence of convexity of the i7-function, and of the property of the map 
(EOl): 
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Now, how to lift general quasi-equilibria (and, consequently, also the Grad approxima- 
tions) to a kinetic model? The answer to this question was given in the Ref. j22I. The 
kinetic model for a quasi-equilibrium approximation /(M) has the form: 



dJ + vVf = --(/ - /(M(/))) + g(/(M(/)), /(M(/))). 



(53) 



Here f{M{f)) is the natural map f f{M), 



I 



mkif-f{Mif)))dv = 0, k = l,...,n, 



(54) 



and thus the first term in the right hand side of equation is just BGK-like, whereas 
the second term, function Q(/(M(/)), /(M(/))) is the true (Boltzmann) collision inte- 
gral, evaluated on the quasi-equilibrium manifold. The latter is crucial: Unlike the true 
Boltzmann collision integral Q{f,f) which can take values in the entire phase space of 
distribution function, Q{f{M{f)), f{M{f))) is allowed to take values only on a relatively 
thin subset known a priori, and can be thus pre-computed to the explicit function of M 
and V (see Ref. |22j for examples). If the quasi-equilibrium f{M) consists only of the local 
Maxwelhans, then Q{f{M{f)), f{M{f))) equals to zero, and we get back the BGK-model. 
In all other cases, the second term in the kinetic model ()53p is essential: If it is omitted in 
equation ()53|) then the zero of the resulting collision integral is the whole quasi-equilibrium 
manifold f{M), and not its local Maxwellian submanifold, unlike the case of the Boltzmann 
collision integral. 



The if-theorem for kinetic models has the following structure 122]: Let us compute 



Function ctbgk is the contribution from the BGK-like term in equation (j^Hj) . and it is 
always non-positive, again due to the property of the map / —>■ f{M) ()54|) . The second 
contribution, ag may be not sign-definite if / is taken far away from the quasi-equilibrium. 
However, one proves j22 that there always exists a non-empty neighborhood of the quasi- 
equilibrium manifold, where < (this is almost obvious: On the quasi-equilibrium 
manifold aQ{f{M)) is the entropy production due to the true Boltzmann collision integral). 
Thus, if the relaxation towards quasi-equilibrium states is fast enough (r is sufficiently close 
to zero), the net entropy production inequality holds, cr = ctbgk + o"q < 0. 

Further simplification of the models (jSHj) are possible. Let us mention here two of them: 



ctrgk + ctq, 




(55) 
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First, instead of function Q{f{M{f)), f{M{f))), we can use 



pQumf)), fmm = E -jkr mmu)), 




k M=M{f) 




(56) 



That is, instead of the true coUision integral Q, we take only its quasi-equilibrium projec- 
tion, PQ. The simplification here is that the velocity dependence is now accumulated only 
in the quasi-equilibrium distribution, and not in the function Q(/(M(/)), /(M(/))): 



Second, and probably the last simplification occurs if one uses the BGK collision integral 
Qbgk with a different relaxation time, say 6, instead of the Boltzmann collision 

integral: 



Here denotes the fc-order moment of the local Maxwellian. 

As a final comment here, the family of the kinetic models reviewed in this section use 
the natural map f f{M) (j54|) of the quasi-equilibrium approximations. Different maps 
/ — > f{M) which do not obey (j54|) were used recently to establish BGK-type models for 
various quasi-equilibrium approximations ^j. 

4 Lattice Boltzmann and other minimal kinetic mod- 
els 

The past decade has witnessed a rapid development of minimal kinetic models for nu- 
merical simulation of complex macroscopic systems. The lattice Boltzmann method is 
particularly valuable minimal extension of the Navier-Stokes equation finding increasingly 
more applications in computational fluid dynamics. Some relation of the lattice Boltzmann 
method to Grad's method was indicated in |56j: Once the Grad method is supplemented 
by the Gauss-Hermite quadrature in the velocity space, the moment system can be rewrit- 
ten in the form of a discrete-velocity model, that is, it becomes amenable to effective 
numerical implementation. Recently, a quasi-equilibrium version of this construction was 
established, in which the quadrature is done not on the distribution function but on the 
entropy functional jH IH] . Quite remarkably, the quasi-equilibrium perspective on the lat- 
tice Boltzmann method results in its refinement known as the entropic lattice Boltzmann 
method 1121 ■ Here we review the entropic lattice Boltzmann method (ELBM) for 
hydrodynamics. 




(57) 




18 



We start with a generic discrete velocity kinetic model. Let fi{x,t) be populations 
of the D-dimensional discrete velocities Vi, i = 1, . . . ,nd, at position x and time t. The 
hydrodynamic fields are the linear functions of the populations, namely 

"d 

v„ v^}f, = {p, pu, pDT + pu^}, (59) 

i=l 

where p is the mass density of the fiuid, pu is the D-dimensional momentum density vector, 
and e = pDT + pu'^ is the energy density. In the case of isothermal simulations, the set 
of independent hydrodynamic fields contains only the mass and momentum densities. It 
is convenient to introduce nd-dimensional population vectors /, and the standard scalar 
product, {f\g) = ^"=i^i2/i- For example, for almost-incompressible hydrodynamics (leav- 
ing out the energy conservation), the locally conserved density and momentum density 
fields are written as 

(l|/)=p, iv^\f)=pUo.. (60) 

Here 1 = Va = {via}^=i, and a = 1, . . . , D, where D is the spatial dimension. 

The construction of the kinetic simulation scheme begins with finding a convex function 
of populations H (entropy function), which satisfies the following condition: If /^^(p, n) 
(local equilibrium) minimizes H subject to the hydrodynamic constraints (equations 
or ()60|) ). then f'^'^ also satisfies certain restrictions on the higher-order moments. For 
example, the equilibrium stress tensor must respect the Galilean invariance, 

rid 

ViaVif3fi'^{p, U) = pcl5af3 + pUaUp. (61) 

i=l 

The corresponding entropy functions for the isothermal and the thermal models were found 
in 1121 El El E] , and are given below (see section E.U.ll and Tabled)). For the time being, 
assume that the convex function H is given. 

The next step is to obtain the set of kinetic equations, 

dtfi + Vio^dafi = Ai. (62) 

Let mi, . . . , m„^ be the nd-dimensional vectors of locally conserved fields. Mi = (mj|/), 
i = l,...,nc. He < fid- The nd-dimensional vector function A (collision integral), must 
satisfy the conditions: 

(mil A) = 0, i = 1, . . . ,nc (local conservation laws), 

(Vi/|A) < (entropy production inequality), 

where Vif is the row- vector of partial derivatives dH/dfi. Moreover, the local equilibrium 
vector /'^^ must be the only zero point of A, that is, A(/'^*^) = 0, and, finally, f^"^ must be 
the only zero point of the local entropy production, (T{f'^'^) = 0. Collision integrals which 
satisfies all these requirements are called admissible. Let us discuss several possibilities of 
constructing admissible collision integrals. 
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BGK model 



Suppose that the entropy function H is known. If, in addition, the local equilibrium 
is also known as an explicit function of the locally conserved variables (or some reliable 
approximation of this function is known), the simplest option is to use the Bhatnagar- 
Gross-Krook (BGK) model. In the case of isothermal hydrodynamics, for example, we 
write 

A = --(/-r^(p(/),n(/))). (63) 
r 

The BGK collision operator is sufficient for many applications. However, it becomes advan- 
tageous only if the local equilibrium is known in a closed form. Unfortunately, often only 
the entropy function is known but not its minimizer. For these cases one should construct 
collision integrals based solely on the knowledge of the entropy function. We present here 
two particular realizations of the collision integral based on the knowledge of the entropy 
function only. 

Quasi- chemical model 

For a generic case of ric locally conserved fields, let g^, s = 1, . . . ,12^ — ric, be a basis of the 
subspace orthogonal (in the standard scalar product) to the vectors of the conservation 
laws. For each vector g^, we define a decomposition gs = gt ~ 9s y where all components 
of vectors gf are nonnegative, and if gf^ 7^ 0, then g^- = 0. Let us consider the collision 
integral of the form: 

A = 5] wsg, {exp {iVH\g:)) - exp {iVH\gt)) } . (64) 

s=l 

Here Wg > 0. By construction, the collision integral (j64|) is admissible. If the entropy 
function is Boltzmann-like, and the components of the vectors g^ are integers, the collision 
integral assumes the familiar Boltzmann-like form. 

Single relaxation time gradient model 

The BGK collision integral (|63|) has the following important property: the linearization of 
the operator (jU^ around the local equilibrium point has a very simple spectrum {0, — 1/r}, 
where is the ric-times degenerate eigenvalue corresponding to the conservation laws, while 
the non-zero eigenvalue corresponds to the rest of the (kinetic) eigenvectors. Nonlinear 
collision operators which have this property of their linearizations at equilibrium are called 
single relaxation time models (SRTM). They play an important role in modelling because 
they allow for the simplest identification of transport coefficients. 

The SRTM, based on the given entropy function H, is constructed as follows (single 
relaxation time gradient model, SRTGM). For the system with Uc local conservation laws. 
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let e^, s = 1, . . . ,nd — Uc, be an orthonormal basis in the kinetic subspace, (miles) = 0, 
and (e^lcp) = 6sp. Then the single relaxation time gradient model is 



nd-ric 

A = — V e,K,pif)iep\VH), (65) 
r ^ — ^ 

s,p=l 

where Kgp are elements of a positive definite (na — ric) x [rid — Uc) matrix K, 

K{f) = C-\f), (66) 
Cspif) = (e,|VVi7(/)|ep). 

Here, VVi/(/) is the rid x "^d matrix of second derivatives, d^H/dfidfj. Linearization of 
the collision integral at equilibrium has the form. 



— e^Cs, (67) 

s = l 



which is obviously single relaxation time. Use of the SRTGM instead of the BGK model 
results in the same hydrodynamics even when the local equilibrium is not known in a closed 
form. Further details of this model and its numerical implementation can be found in Ref. 

It is pertinent to our discussion to explain the term "gradient" appearing in the 
name SRTGM. In Euclidean spaces with the given scalar product, we often identify the 
differential of a function f{x) with its gradient: in the orthogonal coordinate system 
(grad/(a;))j = df{x)/dxi. However, when dealing with a more general setting, one can 
run into problems while making sense out of such a definition. What to do, if there is no 
distinguished scalar product, no preselected orthogonality? 

For a given scalar product ( | ) the gradient grad^/(x) of a function f{x) at a point x 
is such a vector g that {g\y) = D^fiy) for any vector y, where D^f is the differential of 
function / at a point x. The differential of function / is the linear functional that provides 
the best linear approximation near the given point. 

In order to transform a vector into a linear functional one needs a pairing, that means a 
bilinear form (|). This pairing transforms vector g into linear functional {g\: {g\{x) = {g\x). 
Any twice differentiable function f{x) generates a field of pairings: at any point x there 
exists a second differential of /, a quadratic form {Dlf){Ax, Ax). For a convex function 
these forms are positively definite, and we return to the concept of scalar product. Let 
us calculate a gradient of / using this scalar product. In coordinate representation the 
identity (grad/(x) | y)^ = {Dxf){y) (for any vector y) has a form 

T.^^'^^f^-))^^y^ = ll^y^^ (68) 



hence, 



(grad/(x)). = Y.^Dlf)T'^^. (69) 
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As we can see, this grad/(x) is the Newtonian direction, and with this gradient the method 
of steepest descent transforms into the Newton method of optimization. 

Entropy is the concave function, and we define the entropic scalar product through 
negative second differential of entropy. Let us define the gradient of entropy by means of 
this scalar product: (gradj.S'lz)^ = {DxS){z). The entropic gradient system is 

^ = ^(a;)grad^^, (70) 

where ^{x) > is a positive kinetic multiplier. The entropic gradient models fl70|) possesses 
all the required properties (if the entropy Hessian is sufficiently simple). In many cases it is 
simpler than the BGK model, because the gradient model is local in the sense that it uses 
only the entropy function and its derivatives at a current state, and it is not necessary to 
compute the equilibrium (or quasi-equilibrium for quasi-equilibrium models). The entropic 
gradient model has a one-point relaxation spectrum, because near the equilibrium x'^^ the 
gradient vector field ()70|) has an extremely simple linear approximation: d{Ax)/ dt = 
—<f{x'^^)Ax. It corresponds to a well-known fact that the Newton method minimizes a 
positively defined quadratic form in one step. The SRTGM discrete velocity model 

is a particular realization of this construction when the local conservation laws are 
projected out. 



4.0.1 //-functions of minimal kinetic models 

The Boltzmann entropy function written in terms of the one-particle distribution func- 
tion f{x,v) is if = J f In fdv, where v is the continuous velocity. Close to the global 
(reference) equilibrium, this integral can be approximated by using the Gauss-Hermite 
quadrature with the weight 

l^ = (2 7rro)(^/2)exp(-t;V(2ro)). 

Here D is the spatial dimension, Tq is the reference temperature, while the particles mass 
and Boltzmann's constant k-Q are set equal to one. This gives the entropy functions of the 
discrete- velocity models [12113 El, 

//^i;/.ln(A). ,71) 

1=1 ^ ' 

Here, Wi is the weight associated with the i-th discrete velocity Vi (zeroes of the Hermite 
polynomials). The discrete- velocity distribution functions (populations) fi{x) are related 
to the values of the continuous distribution function at the nodes of the quadrature by the 
formula, 

Mx) = Wi{2nToY''/'^ exp{v^/{2To))f{x,v,). 

The entropy functions (f7T|) for various {wi, Vi} are the only input needed for the construc- 
tion of minimal kinetic models. 

With the increase of the order of the Hermite polynomials used in evaluation of the 
quadrature (ffT|l . a better approximation to the hydrodynamics is obtained. The first few 
models of this sequence are represented in Tabled 
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Table 1: Minimal kinetic models j8j. Column 1: Order of Hermite velocity polynomial used 
to evaluate the Gauss-Hermite quadrature; Column 2: Locally conserved (hydrodynamic) 
fields; Column 3: Discrete velocities for D = 1 (zeroes of the corresponding Hermite 
polynomials). For D > 1, discrete velocities are all possible tensor products of the one- 
dimensional velocities in each component direction; Column 4: Weights in the entropy 
formula ()71|). corresponding to the discrete velocities of the Column 3. For D > 1, the 
weights of the discrete velocities are products of corresponding one- dimensional weights; 
Column 5: Macroscopic equations for the fields of Column 2 recovered in the hydrodynamic 
limit of the model. 



1. Order 


2. Fields 


3. Velocities 


4. Weights 


5. Hydrodynamic limit 
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Diffusion 






-VTo 
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Isothermal Navier-Stokes 








1/6 










1/6 




4 


p, pu, e 




l/[4(3 - v^)] 


Thermal Navier-Stokes 






-V3-VqVTo 


l/[4(3 - v^)] 








Vs + VqV% 


l/[4(3 + v^)] 










l/[4(3 + v^)] 





Entropic lattice Boltzmann method 

If the set of discrete velocities forms the links of a Bravais lattice (with possibly several 
sub- lattices) , then the discretization of the discrete velocity kinetic equations in time and 
space is particularly simple, and leads to the entropic lattice Boltzmann scheme. This 
happens in the important case of the isothermal hydrodynamics. The equation of the 
entropic lattice Boltzmann scheme reads 

Mx + c,6t, t + 5t)- f,{x, t) = (5a{f{x, t))A,{f{x, t)), (72) 

where St is the discretization time step, and (3 G [0, 1] is a fixed parameter which matches 
the viscosity coefficient in the long-time large-scale dynamics of the kinetic scheme fl72j) . 
The function a of the population vector defines the maximal over-relaxation of the scheme, 
and is found from the entropy condition, 

H{f{x, t) + aA{f{x, t)) = H{f{x, t)). (73) 

The nontrivial root of this equation is found for populations at each lattice site. Equation 
()73|) ensures the discrete-time if-theorem, and is required in order to stabilize the scheme 
if the relaxation parameter f3 is close to one. The geometrical sense of the discrete-time 
i?-theorem is explained in Fig. |21 We note in passing that the latter limit is of particular 
importance in the applications of the entropic lattice Boltzmann method to hydrodynamics 
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Figure 2: Entropic stabilization of the lattice Boltzmann scheme with over-relaxation. 
Curves represent entropy levels, surrounding the local equilibrium /*^^. The solid curve 
L is the entropy level with the value H{f) = H{f*), where / is the initial, and /* is 
the maximally over- relaxed population / -|- a A defined by equation ()7H|1 . The vector A 
represents the collision integral, the sharp angle between A and the vector — Vif refiects 
the entropy production inequality. The point M is the state of minimum of the entropy 
function H on the line segment between / and /*. The result of the collision update is 
represented by the point The choice of f3 shown corresponds to the over-relaxation: 

H{f{(3)) > H{M) but H{f{f3)) < H{f). The particular case of the BGK colhsion (not 
shown) would be represented by a vector Abgk, pointing from / towards f^"^, in which 
case M = f"^. Figure from Ref. 

because it corresponds to vanishing viscosity, and hence to numerically stable simulations 
of very high Reynolds number fiows. 

Entropic lattice BGK method (ELBGK) 

An important simplification occurs in the case of the isothermal simulations when the 
entropy function is constructed using third-order Hermite polynomials (see Table Q): the 
local equilibrium population vector can be obtained in closed form *S]. This enables the 
simplest entropic scheme - the entropic lattice BGK model - for simulations of isothermal 
hydrodynamics. We present this model in dimensionless lattice units. 

Let D be the spatial dimension. For D = 1, the three discrete velocities are 

« = 0,1}. (74) 

For D > 1, the discrete velocities are tensor products of the discrete velocities of these one- 
dimensional velocities. Thus, we have the 9-velocity model for D = 2 and the 27-velocity 
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model for D = 3. The H function is Boltzmann-like: 

3 



D 



i=l 

The weights Wi are associated with the corresponding discrete velocity Vi. For D = 1, the 
three-dimensional vector of the weights corresponding to the velocities ()74|) is 

For D > 1, the weights are constructed by multiplying the weights associated with each 
component direction. 

The local equilibrium minimizes the if-function i]71\i subject to the fixed density and 
momentum, 

3^ 3^ 

^f^ = ^ fi'^io = pUa, a=l,...,D. (77) 

1=1 i=l 

The explicit solution to this minimization problem reads, 



r = P^. n (2 - yiTa^) f ^"°;^|^ ) . (78) 



a=l 



Note that the exponent. Via, in (f7H|l takes the values ±1, and only, and the speed of 
sound, Cg, in this model is equal to 1/a/3. The factorization of the local equihbrium (fTHj) 
over spatial components is quite remarkable, and resembles the familiar property of the 
local Maxwellians. 

The entropic lattice BGK model for the local equilibrium (j78p reads, 

fiix + Vi6t, t + 6t)- fi{x, t) = -f3aU\{x, t) - ft\p{f{x, t)), u{f{x, t))). (79) 

The parameter f3 is related to the relaxation time r of the BGK model (|63|) by the formula, 

and the value of the over-relaxation parameter a is computed at each lattice site from the 
entropy estimate, 

H{f-aif-r'^im=Hif). (81) 

In the hydrodynamic limit, the model (ffn|) reconstructs the Navier-Stokes equations with 
the viscosity 

/z = pc,V = pc'jt (^-l). (82) 



2(3 



The zero- viscosity limit corresponds to (3 ^ 1. 
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Wall boundary conditions 



The boundary (a solid wall) dR is specified at any point x G dR by the inward unit 
normal e, the wall temperature T^aii, and the wall velocity Wwaii- The simplest boundary 
condition for the minimal kinetic models presented above is obtained upon evaluation of 
the diffusive wall boundary condition for the Boltzmann equation with JE] the help of the 
Gauss- Hermite quadrature. Details can be found in jJl We here write out the final 
expression for the diffusive wall boundary condition: 

fi = ^ fr^ ^TSTT r/i Pwaii, Wwaii , ■ n > . (83) 

<0 K^i' ■ "^Jl/i' (Pwall, lAwall) 

Here is the discrete velocity in the wall reference frame, = Vi — Uwaii- 



Numerical illustrations of the ELBGK 

The Kramers problem |^ is a limiting case of the plane Couette flow, where one of the 
plates is moved to infinity, while keeping a fixed shear rate. The analytical solution for 
the slip- velocity at the wall calculated for the linearized BGK collision model jTHj with the 
simulation of the entropic lattice BGK model are compared in Fig. El This shows that one 
important feature of original Boltzmann equation, the Knudsen number dependent slip at 
the wall is retained in the present model. 

In another numerical experiment, the ELBGK method was tested in the setup of the 
two-dimensional Poiseuille flow. The time evolution of the computed profile as compared to 
the analytical result obtained from the incompressible Navier-Stokes equations is demon- 
strated in Fig. m 

The entropic lattice Boltzmann method upgrades the standard lattice BGK scheme 
[HO] to efficient, accurate and unconditionally stable simulation algorithm for high Reynolds 
number flows IH] . As an illustration, we present the result of comparison of the entropic 
lattice Boltzmann scheme versus the accurate spectral element code in the setup of the 
freely decaying two-dimensional turbulence, see Fig. 

The essential difference between the lattice Boltzmann and the much earlier main body 
research on discrete velocity models pioneered by the seminal work of Broadwell is in 
two points: 

• In the lattice Boltzmann, the effort is done on fixing as much as possible of the 
true (known from continuum theory) Maxwellian dependence of relevant higher-order 
moments on the hydrodynamic moments with as minimum of discrete velocities as 
possible, and 

• The space-time discretization to allow for large time steps (of the order of the mean 
free flight rather than of the collision). This is at variance with most of the numerical 
implementation of discrete velocity models using finite difference ideology. 
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Figure 3: Relative slip at the wall in the simulation of the Kramers problem for shear rate 
a = 0.001, box length L = 32, Voo = a x L = 0.032. Figure from Ref. [S], computed by S. 
Ansumali. 

Outlook: Lattice Boltzmann and microflows 

Gas flows at the micrometer scale constitute a major portion of contemporary fluid dy- 
namics of engineering interest. Because of its relevance to the engineering of micro electro- 
mechanical systems (MEMS), the branch of computational fluid dynamics focused on micro 
scale phenomena is often called "microfluidics" fn\ l^j . Microflows are characterized by 
the Knudsen number, Kn, which is defined as the ratio of the mean free path of molecules 
A and the characteristic scale L of variation of hydrodynamic fields (density, momentum, 
and energy). For typical flows in microdevices, Kn ~ A/L varies from Kn <^ 1 (almost- 
continuum flows) to Kn ~ 1 (weakly rarefied flows). Another characteristic property of 
microflows is that they are highly subsonic, that is, the characteristic flow speed is much 
smaller than the speed of sound. This feature is characterized by the Mach number. 
Ma ~ u/cs, where u is the characteristic flow speed, and Cg is the (isentropic) speed of 
sound. Thus, for microflows. Ma ^1. To be more specific, typical flow velocities are 
about 0.2 m/s, corresponding to Ma ~ 10~^, while values of the Knudsen number range 
between 10~^ < Kn < 10~^. Finally, in the majority of applications, microflows are 
quasi- two- dimensional . 

Theoretical studies of gas flows at finite Knudsen number have begun several decades 
ago in the realm of the Boltzmann kinetic equation. To that end, we mention pioneering 
contributions by Cercignani |TE], Sone jSHI- These studies focused on obtaining either exact 
solutions of the stationary Boltzmann kinetic equation, or other model kinetic equations 
in relatively simple geometries (most often, infinite or semi- infinite rectangular ducts), or 
asymptotic expansions of these solutions. 
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Figure 4: Development of the velocity profile in the Poiseuille fiow. Reduced velocity 
Uyix) = Uy/uy^^^ is shown versus the reduced coordinate across the channel x. Solid 
line: Analytical solution. Different lines correspond to different instants of the reduced 
time, increasing from bottom to top. Symbol: simulation with the ELBGK algorithm. 
Parameters used are: viscosity /i = 5.0015 x 10~^ (/3 = 0.9997), steady state maximal 
velocity = 1.10217 x 10~^. Reynolds number Re = 1157. Figure from Ref. 

computed by S. Ansumali. 

While analytical solutions are important for a qualitative understanding of microfiows, 
and also for the validation of numerical schemes, they certainly do not cover all the needs 
of computational fiuid dynamics of practical interest. At present, two CFD strategies for 
microfiows are well established. 

• Equations of continuous fluid mechanics with slip boundary conditions. The simplest 
semi-phenomenological observation about microfiows is the break down of the no-slip 
boundary condition of fiuid mechanics with increasing Knudsen number. Since mi- 
crofiows are highly subsonic, this leads to the simplest family of models, equations of 
incompressible or compressible fiuid dynamics supplemented by slip velocity bound- 
ary conditions (a review can be found in ^1]). This approach, although widely used 
at the early days of microfiuidics, remains phenomenological. Moreover, it fails to 
predict phenomena such as non-trivial pressure and temperature profiles observed by 
more microscopic approaches. 

• Direct simulation of the Boltzmann kinetic equation. On the other extreme, it is 
possible to resort to a fully microscopic picture of collisions, and to use a molecular 
dynamics approach or a simplified version thereof - the Direct Simulation Monte 
Carlo method of Bird (DSMC) [lOj. DSMC is sometimes heralded as the method 
of choice for simulation of the Boltzmann equation, and it has indeed proven to be 
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Figure 5: Snapshots of the vorticity field in the freely decaying 2D turbulence at t = 0, 
t = 1000, t = 5000, t = 20000 (from left to right). Time measured in the lattice units. Eddy 
turnover time teddy ~ 700. The Reynolds number based on the mean initial kinetic energy 
E and the box- length L equals Re = lV2E/i' — 13134. Upper row: Spectral method; 
Bottom row: Entropic lattice Boltzmann method. Both computations were performed on 
the grid of the same size (512 x 512 grid points). Figure courtesy S. Ansumali. 
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robust in supersonic, highly compressible flows with strong shock waves. However, 
the highly subsonic flows at small to moderate Knudsen number is not a "natural" 
domain for the DSMC simulations where it becomes computationally intensive |54j . 

Since semi-phenomenological computations are not reliable, and the fully microscopic 
treatment is not feasible, the approach to CFD of microflows must rely on reduced models 
of the Boltzmann equation. Two classical routes of reducing the kinetic equations are well 
known, the Chapman-Enskog method and Grad's moment method. The Chapman-Enskog 
method extends the hydrodynamic description (compressible Navier-Stokes equations) to 
flnite Kn in the form of a Taylor series, leading to hydrodynamic equations of increasingly 
higher order in the spatial derivatives (Burnett's hydrodynamics). Grad's method extends 
the hydrodynamic equations to a closed set of equations including higher-order moments 
(fluxes) as independent variables. Both methods are well suited for theoretical studies of 
microflows. In particular, as was already noted by Grad, moment equations are especially 
well suited for low Mach number flows. 

However, applications of Grad's moment equations or of Burnett's hydrodynamics (or 
of existing extensions and generalizations thereof) to CFD of microflows are limited at 
present because of several reasons. The most severe difficulty is in formulating the boundary 
conditions at the reduced level. Although some studies of boundary conditions for moment 
systems were initiated recently [29j, this problem is far from solved. The crucial importance 
of the boundary condition for microflows is actually expected. Indeed, as the rarefaction 
is increasing with Kn, the contribution of the bulk collisions becomes less signiflcant as 
compared to the collisions with the boundaries, and thus the realistic modelling of the 
boundary conditions becomes increasingly important. 

The entropic lattice Boltzmann method seems to be a promising approach to simulations 
of microflows, and is currently an active area of research iSHl IMl HH]- In contrast 
to Grad's method, ELBM is much more compliant with the boundary conditions (see 
above the diffusion wall approximation, which was also rediscovered in [22], where ELBM 
simulations were tested against molecular dynamic simulations with a good agreement). 
Interested reader is directed to two recent papers [71 |Hj where relations between the Grad 
and the lattice Boltzmann constructions are considered in more detail. 

5 Concluding remarks 

The aim of this review was to give a bird-eye picture of the method of moments pioneered 
by Harold Grad a half century ago. Three relatively new issues pertinent to the question 
what physics is beyond Grad's moment approximation and how to obtain this physics were 
discussed in some detail: "Other variables" (triangle entropy method), invariance correc- 
tions, and lifting Grad's equations to a kinetic model. We believe that further development 
of Grad's approach along the lines indicated here will be beneflcial to emerging flelds of 
fluid dynamics, and this review "will be of value for both engineers and mathematicians ... 
who may attempt to turn the invariance condition equation into rigorous mathematics". 
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as was suggested by the referee of this paper. As per mathematical rigor, the situation is 
at least not hopeless for finite-dimensional systems, such as ordinary differential equations 
of chemical kinetics. However, much more work is needed for infinite-dimensional systems 
like the Boltzmann equation where the present level of mathematical achievements in such 
things as existence and uniqueness of solutions does not allow even to start the rigorous 
talking about construction of invariant manifolds. Some mathematical requirements are 
formulated in the recent book ^7\. 
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